#version 420

/*
 * \note:
 * tensors have to be stored in the following way:
 *
 * diag.x = Dxx     Diagonal elements are stored in one vector
 * diag.y = Dyy
 * diag.z = Dzz
 * offdiag.x = Dyz  Off diagonal elements are stored in another vector
 * offdiag.y = Dxz      where all components are stored in that location
 * offdiag.z = Dxy      that is indexed by the coordinate that is not in the tensor index
*/

float sqr( float a )
{
    return a * a;
}

vec3 getEigenvaluesCardano( in vec3 diag, in vec3 offdiag )
{
  const float M_SQRT3 = 1.73205080756887729352744634151;
  float de = offdiag.z * offdiag.x;
  float dd = sqr( offdiag.z );
  float ee = sqr( offdiag.x );
  float ff = sqr( offdiag.y );
  float m  = diag.x + diag.y + diag.z;
  float c1 = diag.x * diag.y + diag.x * diag.z + diag.y * diag.z
             - ( dd + ee + ff );
  float c0 = diag.z * dd + diag.x * ee + diag.y * ff - diag.x * diag.y * diag.z - 2. * offdiag.y * de;
  float p, sqrt_p, q, c, s, phi;
  p = sqr( m ) - 3. * c1;
  q = m * ( p - ( 3. / 2. ) * c1 ) - ( 27. / 2. ) * c0;
  sqrt_p = sqrt( abs( p ) );

  phi = 27. * ( 0.25 * sqr( c1 ) * ( p - c1 ) + c0 * ( q + 27. / 4. * c0 ) );
  phi = ( 1. / 3. ) * atan( sqrt( abs( phi ) ), q );

  c = sqrt_p * cos( phi );
  s = ( 1. / M_SQRT3 ) * sqrt_p * sin( phi );

  vec3 w;
  w[2] = ( 1. / 3. ) * ( m - c );
  w[1] = w[2] + s;
  w[0] = w[2] + c;
  w[2] -= s;
  return w;
}

float getMajorEigenvalue( vec3 diag, vec3 offdiag )
{
  return getEigenvaluesCardano( diag, offdiag )[0];
}


vec3 getEigenvector( vec3 ABC /*diag without eigenalue i*/, vec3 offdiag )
{
  vec3 vec;
  vec.x = ( offdiag.z * offdiag.x - ABC.y * offdiag.y ) * ( offdiag.y * offdiag.x - ABC.z * offdiag.z );
  vec.y = ( offdiag.y * offdiag.x - ABC.z * offdiag.z ) * ( offdiag.y * offdiag.z - ABC.x * offdiag.x );
  vec.z = ( offdiag.z * offdiag.x - ABC.y * offdiag.y ) * ( offdiag.y * offdiag.z - ABC.x * offdiag.x );

  return normalize( vec );
}

float getFA( vec3 evalues )
{
  float sqrt32 = sqrt( 3. / 2. );
  float I1 = evalues.x + evalues.y + evalues.z;
  float mu = I1 / 3.;
  vec3 deriv;
  deriv.x = ( evalues.x - mu );
  deriv.y = ( evalues.y - mu );
  deriv.z = ( evalues.z - mu );
  float FA = sqrt32 * length( deriv ) / length( evalues );
  return FA;
}

float getFA( vec3 diag, vec3 offdiag )
{
  float cross = dot( offdiag, offdiag );
  float j2 = diag.x * diag.y + diag.y * diag.z + diag.z * diag.x - cross;
  float j4 = dot( diag, diag ) + 2. * cross;
  return sqrt( ( j4 - j2 ) / j4 );
}

